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Abstract:  The  numerical  procedure  associated  with  the  inverse  method 
of  Peggion  {1987)  is  described.  The  model  equations  are  formulated  on  a 
spherical-shell  coordinate  system.  The  method  is  applied  to  a  limited  area 
of  the  basin  of  the  Greenland-Iceland-Norwegian  ^GlN^Sea,  using  clima¬ 
tological  and  historical  data  from  the  Generalised  Digital  Environmental 
Model/ of  the  Naval  Oceanographic  Office,  Washington  D.C. 
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1.  Introduction 

An  inverse  method  designed  for  calculation  of  absolute  velocities  from  hydrographic 
data  was  described  recently  as  a  preliminary  study  of  an  investigation  intended  to 
provide  an  understanding  of  the  large-scale,  climatic  features  of  the  ocean  (Peg- 
gion,  1987;  hereafter  referred  to  as  Part  I).  Although  the  method  has  a  large  range 
of  applications,  it  has  been  developed  with  particular  attention  to  the  Greenland 
Iceland- Norwegian  (GIN)  Sea,  and  to  the  aims  of  the  SACLANTCEN  research 
program  in  that  region  We  present  here  the  description  of  the  numerical  proce¬ 
dure,  specifically  designed  to  apply  the  inverse  method  of  Part  1  to  the  GIN  Sea 
area. 

In  order  to  free  the  theoretical  setting  from  complex  mathematical  frameworks, 
the  model  of  Part  I  was  referred  to  a  cartesian  coordinate  system  centred  at  some 
reference  latitude  and  longitude  that  locally  approximates  the  spherical  shape  of  the 
Earth  In  this  metric  th"  effects  of  sphericity  are  retained  by  approximating,  /,  the 
radial  component  of  planetary  vorticity,  with  a  linear  function  of  the  latitudinal 
coordinate.  This  approximation,  known  as  the  , 1-plane  approximation,  is  valid 
for  regions  of  limited  horizontal  extent  (compared  to  the  Earth's  radius)  located 
at  low.  and  mid  latitudes.  At  high  latitudes  or  in  polar  seas  the  approximation 
breaks  down,  and  different  approaches  are  necessary.  For  this  reason,  the  numerical 
procedure  described  in  this  report ,  is  referred  to  a  spherical-shell  coordinate  system. 
Since  the  scale  of  the  vertical  motion  is  small  in  comparison  with  the  radius  of  the 
earth  r,  the  distance  from  the  centre  of  the  earth  is  replaced  by  the  constant  r 
alone,  and  the  vertical  variable  indicates  the  vertical  distance  above  the  surface. 

The  use  of  the  different  coordinate  system  does  not  alter  the  physics  and  validity 
of  the  model  described  in  1’art  1  It  merely  implies  a  different  formal  presentation 
of  the  model  equations,  which  does  not  affect  the  assumptions  and  conclusions  of 
Part  I  Thus  although  in  Sect  2  we  rewrite  the  fundamental  equations  of  Part  I 
in  the  new  coordinate  system,  we  retain  the  same  notations  of  Part  I  as  much  as 
possible  in  order  to  maintain  the  similarities  between  the  two  mathematical  models. 

Section  It  contains  the  executive  summary  of  the  numerical  code.  Section  4  includes 
some  applications  of  the  method  in  the  GIN  Sea  region,  using  climatological  and  his¬ 
torical  data  from  the  Generalized  Digital  Enviromental  Model  (GDEM)  developed 
by  Tr  Davis  and  described  in  Colburn  et  a)  (1980)  Finally,  Serf.  5  summarizes 
anti  discusses  this  study 
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2.  Reformulation  of  the  model  equations 

The  main  equations  of  the  inverse  method  of  Part  I  are  formulated  h^re  in  a 
spherical-shell  coordinate  system  {<£,  A ,:)  where  <p  is  the  latitude,  A  is  the  lon¬ 
gitude,  and  z  is  the  vertical  coordinate  and  is  zero  at  a  bottom  reference  level  and 
increases  upwards.  Thus  the  model  equations  (2.1)  of  Part  I  (hereafter  referred  to 
as  (1-2.1 )  are 


fv 


1 


PqV  cos  <t> 


/“  = - p*. 

PoT 


o  =  -Pz  -  UP, 


- -  (ux  +  (cos^e)*)  4-  U>2  =  0, 

r  cos  <p 


- - upx  f  -vp#  4  wp,  -  0, 

r  cos  <p  r 


f  2f!  sin  </>, 


(2.1a) 


(2.1b) 


(2.1c) 


(2. Id) 

(2.1c) 

(2. If) 


where  r  is  the  radius  of  the  earth,  the  subscripts  <p.  A,  and  2  denote  partial  dif 
ferentation,  ti  -  rcos^dA/df  is  the  eastward  velocity,  r  rd>f>/dt  is  the  poleward 
velocity,  and  w  is  the  vertical  velocity.  The  variable  p  is  the  hydrostatic  pressure 
associated  with  the  density  distribution  p,  and  po  is  a  density  constant  of  reference; 
/  is  the  Coriolis  parameter;  g  is  the  gravitational  acceleration;  and  f!  is  the  angular 
speed  of  rotation  of  the  earth  The  bottom  of  the  ocean  is  taken  to  be  at  2  =  fi(A,0), 
and  the  mean  sea  surface  displacement  at  the  reference  level  2  =  H.  Equal  ions  ( 2. 1 ) 
are  satisfied  in  the  domain  V'  =  {0O  <  <t><  0i.  A0  <  A  <  A,,  h  <  z  <  H}. 


Following  the  formulation  of  Part  I,  the  velocity  field  is  decomposed  as  follows. 

(u, !>,«>)  =  (u'  +  uh,v'  +  t'b.U'1  4-  ll’b) 
where  (u1,  t’’,  w')  are  the  solution  of 


/«■’ 


P  oc 


i _ r 

cos  <t>  Jh 


vdc. 


/»' 


Par  Jh 


L 


Pf  d:, 


fw't  r,  fiv'  -  Y, 


(2.2a) 


2  - 


SACLANTCEN  SM-204 


and  the  bottom  velocities  (tib,  eb,  «’h)  are  given  bv 

i  .  1 


/«  b 


PqT  COS  4> 


(Pbh 


/»b  =  — (pb)*. 
Par 


/(u'b)*  “  /di’b  +  V,  (2.2b) 

where  pb  is  the  pressure  at  the  bottom  associated  with  the  geostrophic  velocites 
"b,  t’b;  P  =  U/r\  and  Y  -  g/(por2  cos  <£)(/>*(  z  =  /i)h*  -  p*(z  -  h)h\).  Here 
the  parameter  /3  is  assumed  not  to  be  constant  as  in  Part  I,  but  as  a  function  of 
latitude. 


From  Eqs.  (2.2),  the  vertical  velocities  «>'  and  trb  are  specified  as  follows: 


-7  /> 

I«’b  =  7  f  - - - »b  l>\  I  -I'/!,-  '  d|C 

f  \  t  cos  <p  r 


(2.3a) 


(2.3b) 


Finally,  the  master  equation  (1-2  8)  for  the  pressure 
follows: 

1  1  , 

- 0  (/>!,  I\  ■  -  I'  (  p|,  )  ;• 

r  cos  o  r 


at  the  bottom  is  specified  as 
h\  (2.4a) 
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(2.4b) 

Following  the  formulation  of  Sect.  3  of  Part  I,  the  variational  principle  (1-3.4)  in 
the  new  coordinate  system  is  transformed  into 

J(p)  =  (  f  f  ( — - — ap\  ^  -hp#  -  r2  cos  <pd</>d-V  d;  )  -  min,  (2.5) 

\A  Jo  \  r  cos  4>  r  /  ) 

where  D  —  {0o  <  <P  <  <t>\ ,  'n  S  '  -\i}-  The  existence  and  uniqueness  of  the 
solution  of  (2.5)  depends  on  the  assignement  of  boundary  conditions.  In  similarity 
with  the  model  formulation  of  Part  I,  we  require  that  the  pressure  at  the  bottom  to 
be  zero  on  dD  Thus,  the  Euler's  equation  associated  with  the  variational  problem 
(2.5)  is  given  by 


l-'(p)  0 

P  -  0 


3 


on  dD, 


(2.6a) 
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Once  more,  the  assigned  boundary  conditions  in  (2. Ha)  are  arbitrary  and  do  not 
restrict  the  validity  of  the  model. 
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3.  The  numerical  program 


The  numerical  scheme  chosen  in  the  treatment  of  the  model  equations  is  baser! 
on  the  entred  difference  method,  using  a  staggered  grid  as  illustrated  in  Tig.  1 
AM;  i;  .  4 hr  assumption  of  geostrophic  balance  suggests  the  use  of  a  grid  where 
th  location  <u  the  u  and  t>  velocities  are  inverted,  we  have  chosen  the  ('-grid  of 
(irammeltvedt  ( 1969)  1  ‘cause  it.  is  the  one  most,  often  used  in  large-scale  dynamics 
models.  Furthermore,  since  measurements  ■  .  seawater  properties  are  mote  frequent 
and  more  reliable  for  the  upper  levels  of  the  ocean,  a  variable  vertical  grid  in  which 
th''  distance  bet  ween  two  levels  is  a  function  of  height  i «  int;  xliiced.  To  preserve  tin* 
accuracy  of  the  -heme,  valuer  of  anv  given  function  at  any  mesh  point  between 
two  consecutive  levels  (not  necessarily  the  middle  point)  is  computed  by  linear 
interpolation. 


P i  .  j+1 


•  •  • 

o._i .i  u..,.,  n,  ,  u,  ,  p,., 


I-  V  |. 

I I- 1. 1-1  1.1-1  1-1-1 

p\ .  j  —  1 


Fig.  I:  Numerical  grid  showing  horizontal  positions  where  density  and 

constituent  distibutions  are  ssumed  to  have  been  measured,  and  the  positions 
where  either  the  bottom  pressure,  p,  and  the  velocities  w,  and  v  are  computed 
The  range  of  the  horizontal  indices  is  as  follows. 
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Assume  that  the  density  and  constituent  distributions  are  known  on  a  horizontally 
regular  mesh  of  (AW  t  1)  x  (A’V  f  1)  -  X  7  points,  where  (AW  f  l).  (A'V  t  1), 
and  NZ  are  the  number  of  stations  on  the  A,  0.  and  r  axes,  respectively.  When 
the  domain  of  interest  does  not  include  land  points  the  observational  grid  contains 
AW  x  A’V  inner  mesh  points  at  which  the  bottom  pressure  might  be  computed. 
Because  of  the  presence  of  cross-derivative  terms,  the  numerical  formulation  of 
the  problem  (2.fi)  involves  nine  adjacent  bottom-pressure  points.  To  preserve  the 
accuracy  of  the  centered  difference  scheme,  the  coefficient  of  the  differential  operator 
(2. da)  are  computed  at  the  grid  points  as  illustrated  in  Fig.  2.  Finally,  when  no  land 
points  are  present,  tin*  variational  principle  (2.5)  is  reduced  to  a  set  of  AW  >  Ar  V 
equations  for  an  equal  number  of  unknowns: 

|  V  A  Xij  I  )  •  m 

P{  /’.r,  ,  l  rn  ] ,  ,  A  A'  (il  l  ! 

1;  I . W. 

The  boundary  conditions  applied  *'*r  the  solution  of  the  differential  equation  (2.h) 
are  generalized  as  follows.  In  tin*  presence  -t  open  boundaries,  values  of  the  pressure 
pt,  are  assigned  In  the  presence  of  rigid  lateral  boundaries,  the  no-slip  boundary 
condition  t  i  .*•  u  r  tit  is  given  fbi*  implies  that  the  normal  rierivative  of 
ph  vanishes  at  the  coasts  Finally,  tie  lifbTnjit ia!  problem  may  be  written  in  the 
traditional  matrix  form: 

M'  H.  (3.2) 

Finding  the  solution  of  a  linear  alg»hrai.-  equation  system  of  large  order  might  be 
a  challenging  numerical  task  While  in  principle  there  are  standard  methods  for 
solving  such  problems,  the  difficulties  are  totally  practical  and  stem  from  the  cost 
of  computer  calculations  and  tin*  loss  of  accuracy  by  runofr  errors.  Thus,  Eq.  (3.2) 
was  solved  using  solver  routines  from  mathematical  computer  libraries  available  at 

S  ACTANT TEN 

The  input  variables  necessary  for  the  application  of  the  program  may  be  categorized 
as  follows 

(i)  Data  input  this  includes  seawater  property  distribution  fields,  obtained  by 
analysis  of  direct  measurements, 

(ii)  l Wjg/i ring  functions-  as  described  by  (l  3.10b)  and  specifier!  as  input  be- 
cause  they  are  directly  related  to  the  definition  of  the  data  distributions; 

(iii)  (isometry  specification  describing  the  basin  geometry  in  terms  of  the  Cori¬ 
olis  parameter  and  topography 

The  constituent  distributions  ami  topograpliv  ar*  sp*-cilm  1  at  the  sanm  column  of 
the  grid  system. 
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Fig  2:  Numerical  k  r  i  *  I  showing  the 

locations  where  the  bottom  pressure  ami 
the  coefficient  of  the  differential  Hq.  (2.6a) 
are  computed.  The  range  of  the  indeces  of 
the  coefficients  is  as  follows: 
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The  final  forcing  of  (2.f>a)  is  also  computed 
at  the  pressure  points 
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4.  Applications 

In  this  section  we  apply  the  method  to  a  portion  of  the  GIN  Sea  region  Although 
one  would  like  to  attempt  a  solution  for  the  global  general  circulation  of  the  whole 
region,  at  the  present  time  it  seems  appropriate  to  deal  with  only  a  limited  amount 
of  data.  Hydrographic  stations  are  taken  from  the  GDEM  climatological  and  his¬ 
torical  data.  More  specifically:  the  salinity  and  temperature  distributions  for  the 
summer  and  winter  seasons  are  considered;  data  are  given  on  0..V'  •  I }  ,ri  horizontal 
grid  at  24  reference  levels,  density  distribution  is  computed  from  the  equation  of 
state  of  Millero  and  Poisson  (1981);  Since  the  mode!  formulation  neglects  bound¬ 
ary  layer  dynamics,  as  was  observed  in  Part  I,  only  data  below  100  m  depth  are 
considered. 

Figure  3  shows  the  bathymetry  of  the  domain  of  the  numerical  experiments.  The 
region  is  crossed  in  its  North-South  extension  by  a  ridge  dividing  the  western  shallow 
Iceland  Plateau  from  the  eastern  deep  Norwegian  Basin.  The  region  is  also  crossed 
by  the  northern  portion  of  the  Polar  Front.  This  front  separates  the  warm  high- 
salinity  waters  of  the  Norwegian  Basin  which  are  of  North  Atlantic  origin  (the 
Norwegian  Atlantic  Water)  from  the  polar  cold  low-salinitv  waters  of  the  Ireland 
Plateau  (the  Polar  Sea  Water).  Because  of  these  features,  the  domain  is  suitable 
for  testing  the  numerical  model  in  the  presence  of  strong  gradients  of  both  the 
topography  and  the  constituent  distributions 

4.1.  THE  SUMMER  SEASON 


Figure  4  shows  the  constituent  distributions  at  the  upper  level  of  the  numerical 
model.  As  expected,  temperature  is  more  sensitive  to  the  mixing  and  diffusion 
processes  than  salinity.  Moving  northwards,  the  temperature  front  loses  strength, 
but  the  salinity  distribution  is  a  more  uniform  function  of  latitude.  The  resultant 
density  distribution  indicates  that  the  Polar  Front  turns  eastwards  at  around  the 
69th  parallel.  A  latitudinal  section  of  the  constituents  (Fig.  5)  shows  that  the 
Norwegian  Atlantic  Water  is  dynamically  constrained  within  the  Norwegian  Basin 
by  the  topographic  features. 

The  velocity,  resulting  from  the  inversion,  follows  this  same  pattern  Where  the 
salinity  and  temperature  gradients  are  both  strong,  the  surface  poleward  jet  is 
more  marked  (Fig.  6a).  However,  in  the  northern  part  of  the  region,  the  jet  turns 
eastwards,  along  the  edge  of  a  meander  of  Polar  Water  The  dynamics  are  even 
less  clear  on  the  eastern  side,  perhaps  due  to  the  fact  that  all  the  constituents 
exhibit  evidence  of  a  gyre  of  Polar  Water,  but  at  different  locations.  Below  the 
therniocline  the  velocity  pattern  is  quite  different  (Fig  6b);  there  is  a  thin  and 
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weak  poleward  jet  following  the  eastern  side  of  tin*  ridge,  hut  tin*  general  How  is 
everywhere  southwards.  A  latitudinal  section  of  the  poleward  velocity  indicates 
that  the  polar  jet  is  mainly  confined  to  the  upper  levels,  with  the  tl«>w  reversal 
occurring  at  about  *100  in  from  the  surface  On  the  Iceland  side,  the  <  inulation  is 
more  barot ropic  with  the  gyre  extending  down  t < >  the  bottom  (fig  *•»•). 


4.2.  THK  VVINTKH  SEASON 

I  luring  the  winter  season,  the  polar  In  >nt  appears  to  be  m<  *re  uniform  with  latitude 
than  during  the  summer  months  (Fig.  7)  AH  the  constituent  distributions  sip. vs 
the  location  of  maximum  gradient  at  around  0  W  Moreover,  a  latitudinal  cross- 
section  indicates  that  the  region  is  horizontally  and  vertically  more  homogeneous 
than  during  the  summer  season.  In  winter,  the  North  Atlantic  Water  is  mlder 
and  less  salty,  smoothing  the  differences  with  the  Polar  Water,  which  is  basirallv 
unaltered  | Fig.  8). 

The  resultant  velocity  field  indicates  that  th“  sumnu  r  poleward  p-t  has  spread  in 
a  wide  current  which  covers  the  whole  eastern  side  «.f  tin*  region,  and  penetrates 
into  the  Ireland  Plateau  in  the  northern  part  Ibpiallv.  tin*  southern  ll«*w  is  un-re 
evident  and  expanded  in  the  western  region  (Fig  '.'at  Below  rhe  t  henno<Tme.  at 
about  1000  m  depth,  the  western  southward  (low  maintain-  its  basic  features  at 
the  surface,  while  the  dynamics  below  the  poleward  current  is  nior*-  complex  In 
the  proximity  of  t lie  ridge,  at  the  southern  boundary,  'b'-p-  is  still  a  weak,  thin 
poleward  current  which  goes  northwards  along  <!i  western  side  of  the  ridge,  but  is 
reflected  back  on  the  eastern  side  and  feed*  tin  return  deep  How  of  the  Norwegian 
Basin  (Fig  0b).  A  latitudinal  section  indicates  no  substantial  vertical  difference* 
between  the  cold  and  warm  seasons  (  Fig.  Or).  The  p.-l.-wap'  current  is  e*s«-tit iall\ 
confined  in  the  upper  !t00  m  of  the  eastern  region. 
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